<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
 <head>
  <title>slasyf.f</title>
 <meta name="generator" content="emacs 21.3.1; htmlfontify 0.20">
<style type="text/css"><!-- 
body { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default   { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default a { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.string   { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.string a { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.comment   { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.comment a { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
 --></style>

 </head>
  <body>

<pre>
      SUBROUTINE <a name="SLASYF.1"></a><a href="slasyf.f.html#SLASYF.1">SLASYF</a>( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  -- LAPACK routine (version 3.1) --
</span><span class="comment">*</span><span class="comment">     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
</span><span class="comment">*</span><span class="comment">     November 2006
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Scalar Arguments ..
</span>      CHARACTER          UPLO
      INTEGER            INFO, KB, LDA, LDW, N, NB
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Array Arguments ..
</span>      INTEGER            IPIV( * )
      REAL               A( LDA, * ), W( LDW, * )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Purpose
</span><span class="comment">*</span><span class="comment">  =======
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  <a name="SLASYF.19"></a><a href="slasyf.f.html#SLASYF.1">SLASYF</a> computes a partial factorization of a real symmetric matrix A
</span><span class="comment">*</span><span class="comment">  using the Bunch-Kaufman diagonal pivoting method. The partial
</span><span class="comment">*</span><span class="comment">  factorization has the form:
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  A  =  ( I  U12 ) ( A11  0  ) (  I    0   )  if UPLO = 'U', or:
</span><span class="comment">*</span><span class="comment">        ( 0  U22 ) (  0   D  ) ( U12' U22' )
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  A  =  ( L11  0 ) (  D   0  ) ( L11' L21' )  if UPLO = 'L'
</span><span class="comment">*</span><span class="comment">        ( L21  I ) (  0  A22 ) (  0    I   )
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  where the order of D is at most NB. The actual order is returned in
</span><span class="comment">*</span><span class="comment">  the argument KB, and is either NB or NB-1, or N if N &lt;= NB.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  <a name="SLASYF.32"></a><a href="slasyf.f.html#SLASYF.1">SLASYF</a> is an auxiliary routine called by <a name="SSYTRF.32"></a><a href="ssytrf.f.html#SSYTRF.1">SSYTRF</a>. It uses blocked code
</span><span class="comment">*</span><span class="comment">  (calling Level 3 BLAS) to update the submatrix A11 (if UPLO = 'U') or
</span><span class="comment">*</span><span class="comment">  A22 (if UPLO = 'L').
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Arguments
</span><span class="comment">*</span><span class="comment">  =========
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  UPLO    (input) CHARACTER*1
</span><span class="comment">*</span><span class="comment">          Specifies whether the upper or lower triangular part of the
</span><span class="comment">*</span><span class="comment">          symmetric matrix A is stored:
</span><span class="comment">*</span><span class="comment">          = 'U':  Upper triangular
</span><span class="comment">*</span><span class="comment">          = 'L':  Lower triangular
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  N       (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The order of the matrix A.  N &gt;= 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  NB      (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The maximum number of columns of the matrix A that should be
</span><span class="comment">*</span><span class="comment">          factored.  NB should be at least 2 to allow for 2-by-2 pivot
</span><span class="comment">*</span><span class="comment">          blocks.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  KB      (output) INTEGER
</span><span class="comment">*</span><span class="comment">          The number of columns of A that were actually factored.
</span><span class="comment">*</span><span class="comment">          KB is either NB-1 or NB, or N if N &lt;= NB.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  A       (input/output) REAL array, dimension (LDA,N)
</span><span class="comment">*</span><span class="comment">          On entry, the symmetric matrix A.  If UPLO = 'U', the leading
</span><span class="comment">*</span><span class="comment">          n-by-n upper triangular part of A contains the upper
</span><span class="comment">*</span><span class="comment">          triangular part of the matrix A, and the strictly lower
</span><span class="comment">*</span><span class="comment">          triangular part of A is not referenced.  If UPLO = 'L', the
</span><span class="comment">*</span><span class="comment">          leading n-by-n lower triangular part of A contains the lower
</span><span class="comment">*</span><span class="comment">          triangular part of the matrix A, and the strictly upper
</span><span class="comment">*</span><span class="comment">          triangular part of A is not referenced.
</span><span class="comment">*</span><span class="comment">          On exit, A contains details of the partial factorization.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LDA     (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The leading dimension of the array A.  LDA &gt;= max(1,N).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  IPIV    (output) INTEGER array, dimension (N)
</span><span class="comment">*</span><span class="comment">          Details of the interchanges and the block structure of D.
</span><span class="comment">*</span><span class="comment">          If UPLO = 'U', only the last KB elements of IPIV are set;
</span><span class="comment">*</span><span class="comment">          if UPLO = 'L', only the first KB elements are set.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">          If IPIV(k) &gt; 0, then rows and columns k and IPIV(k) were
</span><span class="comment">*</span><span class="comment">          interchanged and D(k,k) is a 1-by-1 diagonal block.
</span><span class="comment">*</span><span class="comment">          If UPLO = 'U' and IPIV(k) = IPIV(k-1) &lt; 0, then rows and
</span><span class="comment">*</span><span class="comment">          columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
</span><span class="comment">*</span><span class="comment">          is a 2-by-2 diagonal block.  If UPLO = 'L' and IPIV(k) =
</span><span class="comment">*</span><span class="comment">          IPIV(k+1) &lt; 0, then rows and columns k+1 and -IPIV(k) were
</span><span class="comment">*</span><span class="comment">          interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  W       (workspace) REAL array, dimension (LDW,NB)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LDW     (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The leading dimension of the array W.  LDW &gt;= max(1,N).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  INFO    (output) INTEGER
</span><span class="comment">*</span><span class="comment">          = 0: successful exit
</span><span class="comment">*</span><span class="comment">          &gt; 0: if INFO = k, D(k,k) is exactly zero.  The factorization
</span><span class="comment">*</span><span class="comment">               has been completed, but the block diagonal matrix D is
</span><span class="comment">*</span><span class="comment">               exactly singular.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  =====================================================================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Parameters ..
</span>      REAL               ZERO, ONE
      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
      REAL               EIGHT, SEVTEN
      PARAMETER          ( EIGHT = 8.0E+0, SEVTEN = 17.0E+0 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      INTEGER            IMAX, J, JB, JJ, JMAX, JP, K, KK, KKW, KP,
     $                   KSTEP, KW
      REAL               ABSAKK, ALPHA, COLMAX, D11, D21, D22, R1,
     $                   ROWMAX, T
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      LOGICAL            <a name="LSAME.109"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
      INTEGER            ISAMAX
      EXTERNAL           <a name="LSAME.111"></a><a href="lsame.f.html#LSAME.1">LSAME</a>, ISAMAX
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           SCOPY, SGEMM, SGEMV, SSCAL, SSWAP
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, MAX, MIN, SQRT
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span>      INFO = 0
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Initialize ALPHA for use in choosing pivot block size.
</span><span class="comment">*</span><span class="comment">
</span>      ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
<span class="comment">*</span><span class="comment">
</span>      IF( <a name="LSAME.127"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( UPLO, <span class="string">'U'</span> ) ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Factorize the trailing columns of A using the upper triangle
</span><span class="comment">*</span><span class="comment">        of A and working backwards, and compute the matrix W = U12*D
</span><span class="comment">*</span><span class="comment">        for use in updating A11
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        K is the main loop index, decreasing from N in steps of 1 or 2
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        KW is the column of W which corresponds to column K of A
</span><span class="comment">*</span><span class="comment">
</span>         K = N
   10    CONTINUE
         KW = NB + K - N
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Exit from loop
</span><span class="comment">*</span><span class="comment">
</span>         IF( ( K.LE.N-NB+1 .AND. NB.LT.N ) .OR. K.LT.1 )
     $      GO TO 30
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Copy column K of A to column KW of W and update it
</span><span class="comment">*</span><span class="comment">
</span>         CALL SCOPY( K, A( 1, K ), 1, W( 1, KW ), 1 )
         IF( K.LT.N )
     $      CALL SGEMV( <span class="string">'No transpose'</span>, K, N-K, -ONE, A( 1, K+1 ), LDA,
     $                  W( K, KW+1 ), LDW, ONE, W( 1, KW ), 1 )
<span class="comment">*</span><span class="comment">
</span>         KSTEP = 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Determine rows and columns to be interchanged and whether
</span><span class="comment">*</span><span class="comment">        a 1-by-1 or 2-by-2 pivot block will be used
</span><span class="comment">*</span><span class="comment">
</span>         ABSAKK = ABS( W( K, KW ) )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        IMAX is the row-index of the largest off-diagonal element in
</span><span class="comment">*</span><span class="comment">        column K, and COLMAX is its absolute value
</span><span class="comment">*</span><span class="comment">
</span>         IF( K.GT.1 ) THEN
            IMAX = ISAMAX( K-1, W( 1, KW ), 1 )
            COLMAX = ABS( W( IMAX, KW ) )
         ELSE
            COLMAX = ZERO
         END IF
<span class="comment">*</span><span class="comment">
</span>         IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Column K is zero: set INFO and continue
</span><span class="comment">*</span><span class="comment">
</span>            IF( INFO.EQ.0 )
     $         INFO = K
            KP = K
         ELSE
            IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              no interchange, use 1-by-1 pivot block
</span><span class="comment">*</span><span class="comment">
</span>               KP = K
            ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Copy column IMAX to column KW-1 of W and update it
</span><span class="comment">*</span><span class="comment">
</span>               CALL SCOPY( IMAX, A( 1, IMAX ), 1, W( 1, KW-1 ), 1 )
               CALL SCOPY( K-IMAX, A( IMAX, IMAX+1 ), LDA,
     $                     W( IMAX+1, KW-1 ), 1 )
               IF( K.LT.N )
     $            CALL SGEMV( <span class="string">'No transpose'</span>, K, N-K, -ONE, A( 1, K+1 ),
     $                        LDA, W( IMAX, KW+1 ), LDW, ONE,
     $                        W( 1, KW-1 ), 1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              JMAX is the column-index of the largest off-diagonal
</span><span class="comment">*</span><span class="comment">              element in row IMAX, and ROWMAX is its absolute value
</span><span class="comment">*</span><span class="comment">
</span>               JMAX = IMAX + ISAMAX( K-IMAX, W( IMAX+1, KW-1 ), 1 )
               ROWMAX = ABS( W( JMAX, KW-1 ) )
               IF( IMAX.GT.1 ) THEN
                  JMAX = ISAMAX( IMAX-1, W( 1, KW-1 ), 1 )
                  ROWMAX = MAX( ROWMAX, ABS( W( JMAX, KW-1 ) ) )
               END IF
<span class="comment">*</span><span class="comment">
</span>               IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 no interchange, use 1-by-1 pivot block
</span><span class="comment">*</span><span class="comment">
</span>                  KP = K
               ELSE IF( ABS( W( IMAX, KW-1 ) ).GE.ALPHA*ROWMAX ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 interchange rows and columns K and IMAX, use 1-by-1
</span><span class="comment">*</span><span class="comment">                 pivot block
</span><span class="comment">*</span><span class="comment">
</span>                  KP = IMAX
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 copy column KW-1 of W to column KW
</span><span class="comment">*</span><span class="comment">
</span>                  CALL SCOPY( K, W( 1, KW-1 ), 1, W( 1, KW ), 1 )
               ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 interchange rows and columns K-1 and IMAX, use 2-by-2
</span><span class="comment">*</span><span class="comment">                 pivot block
</span><span class="comment">*</span><span class="comment">
</span>                  KP = IMAX
                  KSTEP = 2
               END IF
            END IF
<span class="comment">*</span><span class="comment">
</span>            KK = K - KSTEP + 1
            KKW = NB + KK - N
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Updated column KP is already stored in column KKW of W
</span><span class="comment">*</span><span class="comment">
</span>            IF( KP.NE.KK ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Copy non-updated column KK to column KP
</span><span class="comment">*</span><span class="comment">
</span>               A( KP, K ) = A( KK, K )
               CALL SCOPY( K-1-KP, A( KP+1, KK ), 1, A( KP, KP+1 ),
     $                     LDA )
               CALL SCOPY( KP, A( 1, KK ), 1, A( 1, KP ), 1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Interchange rows KK and KP in last KK columns of A and W
</span><span class="comment">*</span><span class="comment">
</span>               CALL SSWAP( N-KK+1, A( KK, KK ), LDA, A( KP, KK ), LDA )
               CALL SSWAP( N-KK+1, W( KK, KKW ), LDW, W( KP, KKW ),
     $                     LDW )
            END IF
<span class="comment">*</span><span class="comment">
</span>            IF( KSTEP.EQ.1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              1-by-1 pivot block D(k): column KW of W now holds
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              W(k) = U(k)*D(k)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              where U(k) is the k-th column of U
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Store U(k) in column k of A
</span><span class="comment">*</span><span class="comment">
</span>               CALL SCOPY( K, W( 1, KW ), 1, A( 1, K ), 1 )
               R1 = ONE / A( K, K )
               CALL SSCAL( K-1, R1, A( 1, K ), 1 )
            ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              2-by-2 pivot block D(k): columns KW and KW-1 of W now
</span><span class="comment">*</span><span class="comment">              hold
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              where U(k) and U(k-1) are the k-th and (k-1)-th columns
</span><span class="comment">*</span><span class="comment">              of U
</span><span class="comment">*</span><span class="comment">
</span>               IF( K.GT.2 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Store U(k) and U(k-1) in columns k and k-1 of A
</span><span class="comment">*</span><span class="comment">
</span>                  D21 = W( K-1, KW )
                  D11 = W( K, KW ) / D21
                  D22 = W( K-1, KW-1 ) / D21
                  T = ONE / ( D11*D22-ONE )
                  D21 = T / D21
                  DO 20 J = 1, K - 2
                     A( J, K-1 ) = D21*( D11*W( J, KW-1 )-W( J, KW ) )
                     A( J, K ) = D21*( D22*W( J, KW )-W( J, KW-1 ) )
   20             CONTINUE
               END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Copy D(k) to A
</span><span class="comment">*</span><span class="comment">
</span>               A( K-1, K-1 ) = W( K-1, KW-1 )
               A( K-1, K ) = W( K-1, KW )
               A( K, K ) = W( K, KW )
            END IF
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Store details of the interchanges in IPIV
</span><span class="comment">*</span><span class="comment">
</span>         IF( KSTEP.EQ.1 ) THEN
            IPIV( K ) = KP
         ELSE
            IPIV( K ) = -KP
            IPIV( K-1 ) = -KP
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Decrease K and return to the start of the main loop
</span><span class="comment">*</span><span class="comment">
</span>         K = K - KSTEP
         GO TO 10
<span class="comment">*</span><span class="comment">
</span>   30    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Update the upper triangle of A11 (= A(1:k,1:k)) as
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        A11 := A11 - U12*D*U12' = A11 - U12*W'
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        computing blocks of NB columns at a time
</span><span class="comment">*</span><span class="comment">
</span>         DO 50 J = ( ( K-1 ) / NB )*NB + 1, 1, -NB
            JB = MIN( NB, K-J+1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Update the upper triangle of the diagonal block
</span><span class="comment">*</span><span class="comment">
</span>            DO 40 JJ = J, J + JB - 1
               CALL SGEMV( <span class="string">'No transpose'</span>, JJ-J+1, N-K, -ONE,
     $                     A( J, K+1 ), LDA, W( JJ, KW+1 ), LDW, ONE,
     $                     A( J, JJ ), 1 )
   40       CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Update the rectangular superdiagonal block
</span><span class="comment">*</span><span class="comment">
</span>            CALL SGEMM( <span class="string">'No transpose'</span>, <span class="string">'Transpose'</span>, J-1, JB, N-K, -ONE,
     $                  A( 1, K+1 ), LDA, W( J, KW+1 ), LDW, ONE,
     $                  A( 1, J ), LDA )
   50    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Put U12 in standard form by partially undoing the interchanges
</span><span class="comment">*</span><span class="comment">        in columns k+1:n
</span><span class="comment">*</span><span class="comment">
</span>         J = K + 1
   60    CONTINUE
         JJ = J
         JP = IPIV( J )
         IF( JP.LT.0 ) THEN
            JP = -JP
            J = J + 1
         END IF
         J = J + 1
         IF( JP.NE.JJ .AND. J.LE.N )
     $      CALL SSWAP( N-J+1, A( JP, J ), LDA, A( JJ, J ), LDA )
         IF( J.LE.N )
     $      GO TO 60
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Set KB to the number of columns factorized
</span><span class="comment">*</span><span class="comment">
</span>         KB = N - K
<span class="comment">*</span><span class="comment">
</span>      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Factorize the leading columns of A using the lower triangle
</span><span class="comment">*</span><span class="comment">        of A and working forwards, and compute the matrix W = L21*D
</span><span class="comment">*</span><span class="comment">        for use in updating A22
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        K is the main loop index, increasing from 1 in steps of 1 or 2
</span><span class="comment">*</span><span class="comment">
</span>         K = 1
   70    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Exit from loop
</span><span class="comment">*</span><span class="comment">
</span>         IF( ( K.GE.NB .AND. NB.LT.N ) .OR. K.GT.N )
     $      GO TO 90
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Copy column K of A to column K of W and update it
</span><span class="comment">*</span><span class="comment">
</span>         CALL SCOPY( N-K+1, A( K, K ), 1, W( K, K ), 1 )
         CALL SGEMV( <span class="string">'No transpose'</span>, N-K+1, K-1, -ONE, A( K, 1 ), LDA,
     $               W( K, 1 ), LDW, ONE, W( K, K ), 1 )
<span class="comment">*</span><span class="comment">
</span>         KSTEP = 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Determine rows and columns to be interchanged and whether
</span><span class="comment">*</span><span class="comment">        a 1-by-1 or 2-by-2 pivot block will be used
</span><span class="comment">*</span><span class="comment">
</span>         ABSAKK = ABS( W( K, K ) )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        IMAX is the row-index of the largest off-diagonal element in
</span><span class="comment">*</span><span class="comment">        column K, and COLMAX is its absolute value
</span><span class="comment">*</span><span class="comment">
</span>         IF( K.LT.N ) THEN
            IMAX = K + ISAMAX( N-K, W( K+1, K ), 1 )
            COLMAX = ABS( W( IMAX, K ) )
         ELSE
            COLMAX = ZERO
         END IF
<span class="comment">*</span><span class="comment">
</span>         IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Column K is zero: set INFO and continue
</span><span class="comment">*</span><span class="comment">
</span>            IF( INFO.EQ.0 )
     $         INFO = K
            KP = K
         ELSE
            IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              no interchange, use 1-by-1 pivot block
</span><span class="comment">*</span><span class="comment">
</span>               KP = K
            ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Copy column IMAX to column K+1 of W and update it
</span><span class="comment">*</span><span class="comment">
</span>               CALL SCOPY( IMAX-K, A( IMAX, K ), LDA, W( K, K+1 ), 1 )
               CALL SCOPY( N-IMAX+1, A( IMAX, IMAX ), 1, W( IMAX, K+1 ),
     $                     1 )
               CALL SGEMV( <span class="string">'No transpose'</span>, N-K+1, K-1, -ONE, A( K, 1 ),
     $                     LDA, W( IMAX, 1 ), LDW, ONE, W( K, K+1 ), 1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              JMAX is the column-index of the largest off-diagonal
</span><span class="comment">*</span><span class="comment">              element in row IMAX, and ROWMAX is its absolute value
</span><span class="comment">*</span><span class="comment">
</span>               JMAX = K - 1 + ISAMAX( IMAX-K, W( K, K+1 ), 1 )
               ROWMAX = ABS( W( JMAX, K+1 ) )
               IF( IMAX.LT.N ) THEN
                  JMAX = IMAX + ISAMAX( N-IMAX, W( IMAX+1, K+1 ), 1 )
                  ROWMAX = MAX( ROWMAX, ABS( W( JMAX, K+1 ) ) )
               END IF
<span class="comment">*</span><span class="comment">
</span>               IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 no interchange, use 1-by-1 pivot block
</span><span class="comment">*</span><span class="comment">
</span>                  KP = K
               ELSE IF( ABS( W( IMAX, K+1 ) ).GE.ALPHA*ROWMAX ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 interchange rows and columns K and IMAX, use 1-by-1
</span><span class="comment">*</span><span class="comment">                 pivot block
</span><span class="comment">*</span><span class="comment">
</span>                  KP = IMAX
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 copy column K+1 of W to column K
</span><span class="comment">*</span><span class="comment">
</span>                  CALL SCOPY( N-K+1, W( K, K+1 ), 1, W( K, K ), 1 )
               ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 interchange rows and columns K+1 and IMAX, use 2-by-2
</span><span class="comment">*</span><span class="comment">                 pivot block
</span><span class="comment">*</span><span class="comment">
</span>                  KP = IMAX
                  KSTEP = 2
               END IF
            END IF
<span class="comment">*</span><span class="comment">
</span>            KK = K + KSTEP - 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Updated column KP is already stored in column KK of W
</span><span class="comment">*</span><span class="comment">
</span>            IF( KP.NE.KK ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Copy non-updated column KK to column KP
</span><span class="comment">*</span><span class="comment">
</span>               A( KP, K ) = A( KK, K )
               CALL SCOPY( KP-K-1, A( K+1, KK ), 1, A( KP, K+1 ), LDA )
               CALL SCOPY( N-KP+1, A( KP, KK ), 1, A( KP, KP ), 1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Interchange rows KK and KP in first KK columns of A and W
</span><span class="comment">*</span><span class="comment">
</span>               CALL SSWAP( KK, A( KK, 1 ), LDA, A( KP, 1 ), LDA )
               CALL SSWAP( KK, W( KK, 1 ), LDW, W( KP, 1 ), LDW )
            END IF
<span class="comment">*</span><span class="comment">
</span>            IF( KSTEP.EQ.1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              1-by-1 pivot block D(k): column k of W now holds
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              W(k) = L(k)*D(k)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              where L(k) is the k-th column of L
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Store L(k) in column k of A
</span><span class="comment">*</span><span class="comment">
</span>               CALL SCOPY( N-K+1, W( K, K ), 1, A( K, K ), 1 )
               IF( K.LT.N ) THEN
                  R1 = ONE / A( K, K )
                  CALL SSCAL( N-K, R1, A( K+1, K ), 1 )
               END IF
            ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              2-by-2 pivot block D(k): columns k and k+1 of W now hold
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              where L(k) and L(k+1) are the k-th and (k+1)-th columns
</span><span class="comment">*</span><span class="comment">              of L
</span><span class="comment">*</span><span class="comment">
</span>               IF( K.LT.N-1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Store L(k) and L(k+1) in columns k and k+1 of A
</span><span class="comment">*</span><span class="comment">
</span>                  D21 = W( K+1, K )
                  D11 = W( K+1, K+1 ) / D21
                  D22 = W( K, K ) / D21
                  T = ONE / ( D11*D22-ONE )
                  D21 = T / D21
                  DO 80 J = K + 2, N
                     A( J, K ) = D21*( D11*W( J, K )-W( J, K+1 ) )
                     A( J, K+1 ) = D21*( D22*W( J, K+1 )-W( J, K ) )
   80             CONTINUE
               END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Copy D(k) to A
</span><span class="comment">*</span><span class="comment">
</span>               A( K, K ) = W( K, K )
               A( K+1, K ) = W( K+1, K )
               A( K+1, K+1 ) = W( K+1, K+1 )
            END IF
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Store details of the interchanges in IPIV
</span><span class="comment">*</span><span class="comment">
</span>         IF( KSTEP.EQ.1 ) THEN
            IPIV( K ) = KP
         ELSE
            IPIV( K ) = -KP
            IPIV( K+1 ) = -KP
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Increase K and return to the start of the main loop
</span><span class="comment">*</span><span class="comment">
</span>         K = K + KSTEP
         GO TO 70
<span class="comment">*</span><span class="comment">
</span>   90    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Update the lower triangle of A22 (= A(k:n,k:n)) as
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        A22 := A22 - L21*D*L21' = A22 - L21*W'
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        computing blocks of NB columns at a time
</span><span class="comment">*</span><span class="comment">
</span>         DO 110 J = K, N, NB
            JB = MIN( NB, N-J+1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Update the lower triangle of the diagonal block
</span><span class="comment">*</span><span class="comment">
</span>            DO 100 JJ = J, J + JB - 1
               CALL SGEMV( <span class="string">'No transpose'</span>, J+JB-JJ, K-1, -ONE,
     $                     A( JJ, 1 ), LDA, W( JJ, 1 ), LDW, ONE,
     $                     A( JJ, JJ ), 1 )
  100       CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Update the rectangular subdiagonal block
</span><span class="comment">*</span><span class="comment">
</span>            IF( J+JB.LE.N )
     $         CALL SGEMM( <span class="string">'No transpose'</span>, <span class="string">'Transpose'</span>, N-J-JB+1, JB,
     $                     K-1, -ONE, A( J+JB, 1 ), LDA, W( J, 1 ), LDW,
     $                     ONE, A( J+JB, J ), LDA )
  110    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Put L21 in standard form by partially undoing the interchanges
</span><span class="comment">*</span><span class="comment">        in columns 1:k-1
</span><span class="comment">*</span><span class="comment">
</span>         J = K - 1
  120    CONTINUE
         JJ = J
         JP = IPIV( J )
         IF( JP.LT.0 ) THEN
            JP = -JP
            J = J - 1
         END IF
         J = J - 1
         IF( JP.NE.JJ .AND. J.GE.1 )
     $      CALL SSWAP( J, A( JP, 1 ), LDA, A( JJ, 1 ), LDA )
         IF( J.GE.1 )
     $      GO TO 120
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Set KB to the number of columns factorized
</span><span class="comment">*</span><span class="comment">
</span>         KB = K - 1
<span class="comment">*</span><span class="comment">
</span>      END IF
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="SLASYF.585"></a><a href="slasyf.f.html#SLASYF.1">SLASYF</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>
